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The flux and nuclear composition of ultra-high energy cosmic rays depend on the cosmic distri- 
bution of their sources. Data from cosmic ray observatories are yet inconclusive about their exact 
location or distribution, but provide a measure for the average local density of these emitters. Due 
to the discreteness of the emitters the flux and nuclear composition is expected to show ensemble 
fluctuations on top of the statistical variations, i.e. "cosmic variance". This effect is strongest for 
the most energetic cosmic rays due to the limited propagation distance in the cosmic radiation back- 
ground and is hence a local phenomenon. For the statistical analysis of cosmic ray emission models 
it is important to quantify the possible level of this variance. In this paper we present a completely 
analytic method that describes the variation of the flux and nuclear composition with respect to the 
local source density. We highlight that proposed future space-based observatories with exposures 
of O(10 6 km 2 sryr) will attain sensitivity to observe these spectral fluctuations in the cosmic ray 
energy spectrum at Earth relative to the overall power-law fit. 

PACS numbers: 96.50.S-, 98.70.Sa, 13.85.Tp 



I. INTRODUCTION 

The nuclear composition of ultra-high energy (UHE) 
cosmic rays (CRs) is unknown and stands as one of 
the most challenging questions in particle astrophysics. 
Composition measurements can be made directly up 
to energies of about 100 TeV with space-based experi- 
ments pp. For higher energies, the nuclear composition 
must be inferred from the extensive air shower (EAS) 
produced by the primary cosmic ray when it interacts 
in the upper atmosphere [2]. In order to collect the elu- 
sive events above 10 19 7 eV (which present an integrated 
flux of less than 1 event per km 2 per steradian and per 
century) observatories with large apertures and long ex- 
posure time are needed. 

Today, the leading role is played by ground based fa- 
cilities that cover vast areas with particle detectors over- 
looked by fluorescence telescopes. The largest is the 
Pierre Auger Observatory, with a surface detector ar- 
ray of 1600 water Cherenkov tanks covering 3000 km 2 , 
which accumulates annually about 6 x 10 3 km 2 sryr 
of exposure [3J. The more recently constructed Tele- 
scope Array (TA) covers 700 km 2 with 507 scintillator 
detectors [4], which should accumulate annually about 
1.4 x 10 3 km 2 sryr of exposure. 

In the near future, the JEM-EUSO mission will or- 
bit the Earth on board the International Space Sta- 
tion at an altitude of of about 400 km. Whilst in the 
"nadir" mode, the remote-sensing space instrument (with 
±30° field of view) will monitor an area of approximately 
1.3 x 10 5 km 2 , recording video clips of fast UV flashes by 
sensing the fluorescence light produced through charged 
particle interactions. This innovative pathfinder mission 



will observe approximately 6 x 10 4 km 2 sryr annually [5], 
a factor of 10 above Auger. 

At present, the best indicator of the nuclear compo- 
sition is the atmospheric depth at which the shower de- 
velops its maximum size, A max . For a given shower, the 
position of A max depends on the depth of the first inter- 
action in the atmosphere and the interaction step length 
of the cascade development. Thereby, it depends on the 
scattering cross section of the primary particle with air 
and on features of hadronic interactions at high energies. 
The average shower maximum, (A max ), scales approxi- 
mately as \n(E/A), where E is the energy and A is the 
atomic mass of the primary cosmic ray which generated 
the shower [6]. On average, the shower maximum for 
protons occurs deeper in the atmosphere than that for 
the same energy iron nucleus, (Ag^) > (X^ ax ) . 

Further insight is expected to come from the magni- 
tude of the root-mean-square (RMS) fluctuation of A max . 
If a group of showers is selected within a narrow range 
of energies, then fluctuations about the mean of A max 
(or in a parameter that is closely related to A max such 
as the steepness of the lateral distribution function or 
the spread of muon densities measured at ground) are 
expected to be larger for protons than for iron nuclei. 

Unfortunately, because of the highly indirect method 
of measurement, extracting precise information from 
EASs has proved to be exceedingly difficult. The most 
fundamental problem is that the first generations of par- 
ticles in the cascade are subject to large inherent fluc- 
tuations and consequently this limits the event-by-event 
energy resolution of the experiments. In addition, the 
center-of-mass energy of the first few cascade steps is 
well beyond any reached in collider experiments. There- 
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fore, one needs to rely on hadronic interaction models 
that attempt to extrapolate, using different mixtures of 
theory and phenomenology, our understanding of particle 
physics. 

HiRes has presented evidence that the CR composition 
remains proton- like up to the highest energies [7J. The 
TA data are also compatible with proton primaries, but 
the statistics are still limited [8]. In contrast, Auger data 
on the depth of shower maximum, its RMS fluctuation, 
and muon rates at ground level indicate a transition of 
UHE CRs within the energy range 10 18 eV to 10 19 6 eV 
from a light (presumably proton-dominated) towards a 
heavier composition [S]. The apparent contradictory re- 
sults reported by the HiRes, the Telescope Array, and the 
Pierre Auger collaborations suggest that known [TrjlfT5] 
and/or unknown systematic uncertainties may still affect 
the interpretation of EAS data. 

Various features in the CR spectrum can also pro- 
vide indirect evidence for the nuclear composition of 
UHE CRs. The ankle - a hardening of the spectrum 
at I0 185 eV - could be formed naturally by the superpo- 
sition of two power-law fluxes and serves as a candidate 
of the transition between galactic heavy nuclei and extra- 
galactic cosmic ray protons [131 HI] • It nas a ls° been ad- 
vocated that this feature could be well reproduced by a 
proton-dominated power-law spectrum, where the ankle 
is formed as a dip in the spectrum from the energy loss of 
protons via Bethe-Heitler pair production |15U16| . In this 
case extra-galactic protons would already have started to 
dominate the spectrum beyond the 2nd knee, a feature 
which corresponds to a slight softening of the spectrum 
at 10 17 - 7 eV. 

Proton-dominance beyond the ankle is ultimately lim- 
ited by the onset of photopion production on the cosmic 
microwave background (CMB), whereas dominance of a 
heavy composition is restricted by nucleus photodisinte- 
gration through the giant dipole resonance (GDR) [T71 
ITS] - the so called GZK- suppression at around 10 19 7 eV. 
Indeed, a flux suppression in this energy region has been 
observed in the HiRes and Auger data with a high sta- 
tistical significance [TTjH5T] . As noted elsewhere [55J 123] , 
secondary neutrinos and 7-rays of these hadronic inter- 
actions can serve as additional discriminators between 
various CR models. We will return to this interesting 
possibility in our conclusions. 

In this paper we elaborate on the question as to what 
extent the spectral information in the GZK region can be 
used to discriminate between different CR source compo- 
sition models. Due to the strength of the GZK mecha- 
nism the spectrum in this region is dominated by (and 
requires the presence of) local sources [53]. In this case 
the flux from a few CR sources can significantly fluctu- 
ate from a homogeneous distribution that is typically as- 
sumed in CR flux predictions [25] . In contrast to Poisson 
fluctuations in the GZK region [55] the manifestations of 
ensemble fluctuations persist in the limit of large event 



statistics. We will quantify these stochastic fluctuations 
in the following utilizing an analytic solution to the flux 
of CR nuclei derived in Refs. [571 HE]. 

We will start in Sec. [H]with a brief review of the prop- 
agation of CR nuclei and the calculation of the mean 
observed fluxes. We discuss the analytic solution of the 
equation describing the flux from local CR sources in 
Sec. IIIH In Sec. IIVI these results will be used to derive an 
analytic approximation of the flux and mean mass vari- 
ations due to the distribution of sources. We summarize 
our findings in Sec. [V] 



II. PROPAGATION OF COSMIC RAY NUCLEI 

Since cosmic rays are subject to deflections in galac- 
tic and inter-galactic magnetic fields the observed CR 
events do not point directly back to their sources. The 
identification of the CR sources is hence experimentally 
challenging and has so far proved inconclusive. There is 
a general consensus that the sources which are responsi- 
ble of the UHE CR spectrum are of extra-galactic origin. 
These sources are expected to follow a spatially homo- 
geneous distribution and the mean (ensemble- averaged) 
flux of UHE CRs (of type i) follows a set of (Boltzmann) 
continuity equations of the form: 

% = dsiHEYi) + d E {b t Y t ) - rf 4 Y z 

+ dEj^Yj + HQi, (1) 
3 J 

together with the Friedman-Lemaitre equations describ- 
ing the cosmic expansion rate H(z) as a function of 
red-shift z. We follow the usual cosmological concor- 
dance model dominated by a cosmological constant with 
57a ~ 0.7 and a (cold) matter component, Q m ~ 0.3 
where H 2 (z) = Hq [fi m (l + z) 3 + normalized to its 
value today of H ~ 70 kms^Mpc" 1 [50]. The time- 
dependence of the red-shift can be expressed via dz — 
—dt (I + z)H . The first and second terms on the r.h.s. of 
Eq. ([!]) describe, respectively, red-shift and other contin- 
uous energy losses (CEL) with rate b = dE/dt. The third 
and fourth terms describe more general interactions in- 
volving particle losses (i — > anything) with total interac- 
tion rate r* ot , and particle generation of the form j — > i 
with differential interaction rate 7^. The last term on 
the r.h.s. corresponds to the emission rate of CRs of type 
i per co-moving volume, depending on the emission rate 
Qi per source and their density H. 

The two main reactions of UHE CR nuclei during their 
cosmic evolution are photodisintegration 30-37] and 
Bethe-Heitler pair production |38j with the cosmic ra- 
diation background. In addition to the dominant contri- 
bution of the CMB we also include the infra-red/optical 
background from Ref. |39j in our calculation of interac- 
tion and energy loss rates. Photodisintegration is domi- 
nated by the GDR with main branches A —> (A — 1) + N 
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and A —> (A — 2) + 2N where N indicates a proton or 
neutron [31]. The GDR peak in the rest frame of the 
nucleus lies at at about 20 MeV for one-nucleon emis- 
sion, corresponding to -Eq DR ~ A x 2 x e m J v x 10 19 eV 
in the cosmic frame with photon energies e = e me v nieV. 
At energies below 10 MeV there exist typically a num- 
ber of discrete excitation levels that can become signif- 
icant for low mass nuclei. Above 30 MeV, where the 
photon wavelength becomes comparable or smaller than 
the size of the nucleus, the photon interacts via substruc- 
tures of the nucleus. Out of these the interaction with 
quasi-deuterons is typically most dominant and forms a 
plateau of the cross section up to the photopion produc- 
tion threshold at ~ 145 MeV. Bcthc-Hcitlcr pair produc- 
tion can be treated as a continuous energy loss process 
with rate bA(z,E) — Z 2 b p (z, E/A), where b p is the en- 
ergy loss rate of protons [38]. The (differential) photo- 
disintegration rate Ta-*b(E) (ja^b(E, E')) is discussed 
in more detail in Ref. [2"B] . 

The evolution of the spectra proceeds very rapidly on 
cosmic time scales and the diffuse flux of secondary nu- 
clei, J, looks generally quite different from the initial 
source injection spectrum [35] ■ The reaction network of 
nuclei depend in general on a large number of stable or 
long-lived isotopes. If the life-time of an isotope is much 
shorter than its photodisintegration rate it can be effec- 
tively replaced by its long-lived decay products in the 
network ([!]). Typically, neutron-rich isotopes /3-decay to 
a stable or long-lived nucleus with the same mass num- 
ber. In most cases there is only one stable nucleus per 
mass number below 56 Fe with the exception of the pairs 
54 Cr/ 54 Fe, 46 Ca/ 46 Ti, 40 Ar/ 40 Ca and 36 S/ 36 Ar. We fol- 
low here the approach of Ref. [31] (PSB) and consider 
only a single nucleus per mass number A in the decay 
chain of primary iron 56 Fe. This PSB-chain of nuclei 
linked by one-nucleon losses is discussed in more detail 
in Ref. [28]. 

Note, that the set of Boltzmann Eqs. ([l} does not take 
into account the deflection of charged CR nuclei dur- 
ing their propagation through magnetic fields. Magnetic 
scattering can be viewed as a diffusion process that de- 
pends on the particle's Larmor radius J?l = E/{ZeB) ~ 
l.lMpci?Ecv/(Zi?nG) and the coherence length of the 
magnetic field. At energies where the diffusion length be- 
comes larger than the distance to the nearest CR source 
the CR spectrum is expected to be suppressed. It has 
been speculated that for particularly strong inter-galactic 
magnetic fields of strength ~ 1 nG and coherence length 
of ~ 1 Mpc (see e.g. Ref. [30] ), the diffusive propagation 
of CR protons can start to affect the spectrum below 
about 10 18 eV [41] , For heavier nuclei diffusive propa- 
gation can in principle remain important up to higher 
energies due to the dependence i?L oc 1/Z [32] • The re- 
sults of this paper assume that the contribution of inter- 
galactic or galactic magnetic fields can be neglected for 
the calculation of the UHE CR spectrum. 

As an example we show in Fig.[l]the solution of Eq. ([!]) 
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FIG. 1: The spectra from pure iron sources with exponen- 
tial cutoff i?max = 10 12 GeV and power index 7 = 2. The 
green line shows the result of the analytic approximation for 
a homogeneous distribution of sources r m i n < r < r max with 
r m in = 10 Mpc and r max = 1 Gpc. The solid black line is the 
corresponding numerical calculation including redshift scal- 
ing of interaction rates and energies. The dashed black line 
shows a model including sources up to z max = 2. Also shown 
is recent data from HiRes |19| . Auger [5T] and the Telescope 
Array (TA) [43]. 

for an iron source model using a power-law pectral emis- 
sion rate QFe oc E~ J exp(— E / E majyi ) with 7 = 2 and 
B max = 10 21 eV. This model is motivated by a previous 
study [53] and reproduces the Auger data above the ankle 
within systematic uncertainties. The dashed black line 
corresponds to source contributions above r min — 10 Mpc 
extending up to redshift z max = 2 where no red-shift evo- 
lution of the emission is assumed, i.e. % oc 1. The solid 
black line marks the local contribution up to 1 Gpc cal- 
culated with the same method. This local contribution 
where red-shift scaling of energies and interaction rates is 
suppressed can be approximated by an analytic solution 
which is shown as the green solid line in the plot. We 
will discuss this method in the next section. 



III. ANALYTIC SOLUTION 

Photodisintegration and photopion interactions of nu- 
clei happen on timescales much shorter than the Hubble 
scale. Since we are interested in the variation of the av- 
erage flux and mass composition from the distribution of 
local CR sources we will neglect redshift scalings in the 
following. In this case the Green's function of the Boltz- 
mann equations (fTl) can be expressed in a simple analytic 
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form, as discussed in Refs. [2"7l 128) . 

The secondary nuclei produced via photodisintegration 
carry approximately the same Lorentz factor as the ini- 
tial nucleus and the differential interaction rate in the 
set of Eqs. Q can be approximated as ja-^b(E,E') ~ 
T A ^ B (E)S(W - (B/A)E). It is hence convenient to 
express the energy of a nucleus with mass number A 
and red-shift z as AE where E denotes the energy 
per nucleon. Introducing the binned CR flux Fa.i = 
AEiAdF ^(AEi) / 'dE , and corresponding emission rates, 
Qaa = AAEiQA(AEi) we can re-write Eq. Q in the 
compact form 1 



jd r (r 2 F A , 

B<A 



I.AA 



(A,i)-t(B,i)FA,i + ^ ^(B,i)^(A,i)Fb,i 
B>A 

^CEL t^CEL , 



+ r^.i+i-^.i+i - T A l F A ,i , (2) 
where we define the rates: 

pCel _ b A (AE t ) 



A,i 



F(A,i)^(B,i) = r 



AAEi ' 

A^B(AEi) 



(3) 
(4) 



In Ref. [2H] it was shown that the general solution of 
Eq. ([2]), describing the flux from a source at distance r 
can be written in the form 



Airr 2 



lB,j 



c fc=l 

with dimensionless amplitudes 

n c — 1 n c 

A k ( C ) = n r Q ^ Ci+1 / n 1 r 



tot ptot 

Q->c i + i J J_ J_ ^ Cp 1 Cfc 

1=1 P=l(jtk) 



(5) 



(6) 



where we sum over all possible production chains c — 
(ci, . . . , c„ c ) with intermediate nuclei of mass number C 
in the energy bin k - denoted by the doublet c, — (C, k) 
- and endpoints ci = (B,j). The partial width T Cl _,. Cl+1 
includes nucleon-disintegration (T( A .i)-y(B .i)) as we ll as 
CEL (T( A ,i)^(A,i-i))- 

In the following we will apply two more approximations 
to the solution Q to reduce the number of terms. Firstly, 
we will reduce the CEL (last two terms on the r.h.s. of 
Eq. |2|) to the effective loss term 



^CEL 

■ A,i+ 



iFa,i+i - ^a^Faa -> 



b A (AEj) 
AE, 



Fa, 



(7) 



1 This form of the differential equation holds for nuclei heavier than 
beryllium. We can easily compensate for the process 9 Be — ¥ 4 He 
+ 4 He + n of the PSB chain (see Ref. [28]) by re-defining F' A . = 
F Al i/2 for A = 2, 3, 4 and F' A ■ = Fa : % for other nuclei. Similarly, 
the chains with nucleons as a final particle are re-weighted by 
the corresponding multiplicity in the case of intermediate nuclei 
lighter than 9 Be. 



This introduces a relative error of the order 
dE{bEF)/(bF), which is small in the relevant en- 
ergy region of 10 18 — 10 19 eV. And, secondly, we will 
only consider the one-nucleon loss chain of the PSB 
approximation. It was shown in Refs. |27l I28| that this 
is a good approximation for CR nuclei with a large mass 
number, which is the focus of this analysis. 



IV. ENSEMBLE AVERAGE AND VARIATION 



We now discuss the statistical variation of the CR flux 
and composition following Ref. [44]. Herein we assume 
that the probability distribution function (PDF) for local 
(t/Hq <C 1) sources is flat in Euclidean space. Let n s 
sources be distributed between redshift r m - ln and r max . 
The number of sources can then be expressed via the 



(local) source density Hq as n s = "Ho(47r/3)(r 
The PDF of a single source is given by 



max mm- 1 



p{r) 



Ho 



47rr 2 6(r - 



x)e(r n 



-) 



(8) 



The ensemble-average of a quantity A(r%, . . . , r„ 3 ) de- 
pending on the distance of the n s sources can then be 
expressed as 



(A) 



dri 



dr n j>(n) ■ ... -p{r„ s )A. (9) 



The ensemble-average of the local flux of CR nuclei given 
by En s F A ,i(r s ) and Eq. dsb is simply 



A.i] 



'max 

Ho J dr'4ir' 2 F A ,( r '), 



(10) 



Using the abbreviation ^ we can then write this in an 
analytic form as 



c k=l 



(^A, i )=EE^( r c° t )^„ c , 

where we define 



Ho 1 _ T 



(ii) 



(12) 



From the experimental point of view the interesting 
quantities are the ensemble-averaged total flux of nuclei 
iVtot and mean mass number A av at the highest CR en- 
ergies E. The mean total flux is simply given by Eq. ( II ) 
as 



(N tot (E)) = ^2(N A (E)) 



(13) 



We now return to the example of an iron source model 
shown in Fig. [l] The solid green line in Fig. [T] shows 
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Fe source, y= 3. r rajl , = 10 Mpc. r mM = 1 Gpc, = 10" 5 Mpc" : 









3 



























Fc source, y 2. r lllm 100 Mpc. r lu;ix I Gpc, ./4 10 ' Mpc ■ 



£ C R [GeV] 

Fe source, /= 3, i- mLn = 10 Mpc, r max = 1 Gpc, ^ = 10 -5 Mpc -3 









'-111 

rl 


! y 7 " 








flfl 
■ 
11 










1 




(' 


,v> 






1 ^ 





5 


I 




3 






















VPS) 
(JV„) 









£ C R [GeVI 

Fe source, 7= 2, r llJn = 100 Mpc, r ma , = 1 Gpc, J^ = 10~ 5 Mpc" 3 




FIG. 2: The local relative error of the flux (Eq. (161; upper plots in red) and average mass composition (Eq. (18 1; lower plots 
in blue) for a distribution of iron sources with the model parameters indicated above the plots. The green line in the left 
plots indicate the corresponding relative error for the model shown in Fig. [I] All calculations assume a local source density of 
Ho = 10" 5 Mpc" 3 and scale as Hq 1/2 ■ 



the result of the analytic approximation (11). This ap- alytic form as 



proximation agrees with the numerical result (solid black 
line) within a factor two or better depending on energy. 
The relative difference is expected from the onset of red- 
shift scaling for propagation distances of the order of Gpc 
(z ~ 1/4). At low energies this introduces a relative 
upward shift of the flux of ~ 30% for a source model 
with 7 = 2 and H oc 1. This agrees well with the re- 
sult of the calculations. In addition, threshold effects 
that lead to breaks in the spectrum scale with redshift as 
E t h oc 1/(1 + z) 2 and are shifted to lower energies by up 
to ~ 50%. This effect can also be noticed by the relative 
position of the break in Fig. [I] 



(SN Ati SN B j) = (N A ,iN Btj ) - (N Ati )(N Bd ) 

c,c k=l fc=l 

- ~(N A ,i)(N Bd ) , (14) 



where £ is defined as the expression 

7"—. 

,-rT 



C(r) = H J 



dr 



(15) 



This example illustrates the limitations of the analytic 
solution for the calculation of large scale (early-time) con- 
tributions to the CR flux. However, the analytic approx- 
imation provides a convenient description of the nuclei 
cascades that happen on small scales and depend on the 
local source distribution. In particular, it enables us to 
study ensemble-variations of the flux. Using Eq. ([6| we 
can write the variation of the CR flux in an explicit an- 



Notc that the last term in Eq. ( 14 1 is sometimes omitted 



since the number of sources n s is expected to be large, but 
we keep it in our calculations. Based on these definitions 
we can express the relative variation of the total flux via 



the two-point density perturbations ( 14|) as 

f A (E)6N B ( 



2 = v (6N A (E)SN B (E)) 

CT loc 2^ /AT. . 
A.B 



(16) 
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0.14 
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0.01) 


2.5 


173 


0.19 



TABLE I: Estimated source number and adiabatic variation 



defined in Eqs. (19 1 and (21 1. We show results for a power- 
law cosmological evolution as H = Ho(l + z) n with z m i n < 
z < Zmax and for the star formation rate (SFR) according to 
Rcf. 45, 46 (with cutoff z m j n = 0.01). In all cases we assume 
a local density of Ho = 10" 5 Mpc -3 . 



With Eq. ( 13 ) we can also define the mean mass num- 
ber as 



(A av (E)) = 



A(N A {E)) 
(Ntot(E)) ■ 



(17) 



Note, that Eq. ( 17 ) is in the strict sense not the ensemble- 
average but serves as a first order estimator. For small 
fluctuations around the mean value we can approximate 
the relative variation of the mean mass number (17) via 



the two-point correlation function ( 14 ) as 



E 

B,C 



D 



C 



(SN B (E)SN C (E)) 
(Ntot(E))* 



(18) 



a local distribution between 10 Mpc and 1 Gpc. The 
relative errors are significantly reduced as we increase 
the distance to the closest source to 100 Mpc as shown 
in the plots of the last column in Fig. [2] However, this 
is only marginally reproducing the UHE CR spectrum as 
pointed out in Ref. [24]. Note that the green lines in the 
plots of the left column mark the contribution for the 
iron source model with E max = 10 21 eV considered in 
Fig. 

One can also notice from the lower plots of Fig. [2] that 
the relative error of the average mass composition is be- 
low 1% for CR energies below 10 19 5 eV. This energy 
marks the end of the energy region where CR event statis- 
tics allow an inference of the mass composition from CR 
data as described in the introduction. Note that all the 
plots in Fig. [2] show the case of a local source density of 
Ho = 10 -5 Mpc" 3 and levels increase as Hq 1 ^ 2 . Hence, 
even for a density Ho = 10 -6 Mpc -3 , still marginally 
consistent with the data, the ensemble fluctuation on the 
average mass composition probed by present generation 
experiments may be safely neglected. 

So far we have considered the case that T/H <C 1 valid 
at the highest CR energies where we can neglect cosmo- 
logical contributions of the sources and treat the problem 
as effectively local. In the opposite case, T/Hq ^ 1, it is 
also possible to give an analytic expression of the statisti- 
cal variation of the CR flux. We assume that n s sources 
are isotropically distributed between redshift z m - m and 
z m ax with comoving density H(z), and local density Ho- 
The number of sources is then given by 



dz 



4ird 2 c (z)H(z) (19) 



where the comoving volume is Vc{z) = (47r '/3i)d$j(z) 
with comoving distance (in a flat universe) dc(z) = 
Jq dz' /H(z'). The PDF of a single source as a function 
of redshift is then a simple generalization of Eq. (l8|) , 



In Figure [2] we show contour plots of the relative er- 
ror of the flux (red; upper plots) and mass composition 
(blue; lower plots) for the case of iron sources with differ- 
ent model parameters indicated above the plots. The axis 
show the observed CR energy vs. the exponential cutoff 
energy E max . For the calculation we introduced logarith- 
mic energy bins of Alog 10 E/GeV = 0.02 and smoothed 
the result with a Gaussian kernel to account for an ex- 
perimental resolution of 0.1. This procedure smoothes 
out features in the relative variance that are beyond the 
experimental resolution and result from the rapid nuclei 
transitions in the GZK region in combination with the 
relative CR energy shift with the mass of the daughter 
nuclei. 

The results do not strongly depend on the spectral 
index as can be seen for the cases 7 = 2 and 7 = 3 
with otherwise equal parameters that are shown in the 
plots of the first two columns. This set of plots assumes 



p(z) 



1 H{z) 
H(z) n s 



4*cP c (z). 



(20) 



As long as only adiabatic redshift scaling is involved the 
flux of a single source is given by dF/dE = Q((l + 
z)E)/{^-Kd 2 J ). Assuming Q(E) cx E~ J we then obtain 
a relative variation of 



'adi 



jdz P (z) [{i + z)-ydi{z)Y 
[jdz P (z)(i + z)-ydUz)Y 



i 



(21) 



We show values for n s and (c 2 ^) 1 / 2 in Tab.jl] Obviously, 
the statistical variation of the average mass number is 
negligible in this regime. 

The two limiting behaviors of the relative flux variation 



(16) and (21) motivates the following treatment of the 
overall relative flux variation. Using the numerical result 
of the set of Boltzmann Eqs. ([I]) we can calculate the 
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FIG. 3: Left panel: The example of Fig. [T] including the approximate variation of the flux assuming a local source distribution 
Ho — 10~ 5 Mpc~ 3 (dark gray band) and "Ho = 10 -6 Mpc -3 (light gray band). Right panel: For comparison, a proton model 
with similar source parameters. 



average flux for the local (r < 1 Gpc), and global source 
distribution. This corresponds to the solid and dashed 
lines, respectively, in Fig. [T] The flux variation can then 
be approximated via the superposition 



'A' 



(E) ~ (1 - xYai di + 



i i 
x a 



loc 



(E). 



(22) 



where a? is calculated from the analytic approxima- 



tion (16), tr 2 di the adiabatic limit (21) for r min ~ 1 Gpc 
(z min ~ 0.25) and x = (N loc {E)) /(N glohal {E)) from the 
numerical evaluation. 

In the left plot of Fig. [3] we show the result of this 
procedure for the previous example of iron sources shown 
in Fig. fTl The solid line (corresponding to the dashed line 
in Fig.TTl) shows the contribution from the cosmological 
distribution of iron sources. The shaded areas show the 



flux within its variation based on Eq. (22) for a local 
source density Ho = 10~ 5 Mpc -3 (dark gray) and Ho — 
10~ 6 Mpc~ 3 (light gray). A local source density of Ho = 
10~ 5 Mpc -3 is consistent with the absence of "repeaters" 
in CR data [47l [48] . Note, that the relative size of the 

error bands increases as Hq 1 ^ 2 . Values as low as Ho = 
10~ 6 Mpc -3 are still marginally consistent with auto- 
correlation studies of UHE CR nuclei [22] . 

The right plot of Fig. [3] shows the corresponding result 
of a proton source model with similar source parameters. 
The amplitude of the ensemble fluctuation is compara- 
ble to the case of iron and does not serve as a direct 
measure of the source composition. However, the spec- 
tral feature of the GZK-suppression is significantly dif- 
ferent to the case of the iron model. In particular, the 
iron source model experiences a much steeper suppres- 
sion at the upper end of the UHE CR spectrum. This 



poses a challenge to account for the most extreme UHE 
events like the 3 x 10 20 eV Fly's Eye event 50J, if sources 
are too distant [23] . Candidate nearby sources of heavy 
nuclei include starburst galaxies [SI] and ultra-fast spin- 
ning newly-born pulsars [52] [53] . Our method provides 
a tool to distinguish ensemble fluctuations from spectral 
features of the the average source contribution on a sta- 
tistical basis. 

Proposed future space-based CR observatories can 
reach integrated exposures of O(10 6 km 2 sryr) [54). This 
is a factor of about 100 larger than the integrated expo- 
sure reached by present ground-based air shower arrays. 
The statistical error of UHE CR measurements can thus 
be improved by a factor 10. With such a resolution the 
observed UHE CR spectrum can show a significant de- 
viation from the ensemble mean due to the discreteness 
of close-by source contributions. For instance, the spec- 
trum could exhibit "spectral wiggles" beyond the GZK 
suppression. Our method describes a way to quantify the 
bin-to-bin amplitude of these modulations for the case of 
a ultra-high energy cosmic ray nuclei. 



V. CONCLUSIONS 



In this paper we have studied the ensemble fluctua- 
tions of the mean flux and average mass number of UHE 
CR nuclei from the distribution of sources. We have de- 
rived an analytic expression for the relative errors which 
applies to the CR data at the highest energies dominated 
by the local sub-Gpc (low redshift) source distribution. 



For lower energies another analytic approximation can 
be derived through the use of adiabatic scaling of the 
source contributions. This method can be easily general- 
ized to the case of ensemble fluctuations due to the source 
emission parameters such as the spectral index and the 
maximal energy. 

As an illustration, we applied these results to a fidu- 
cial iron source model for which a homogeneous distribu- 
tion of sources had previously been found to successfully 
reproduce the CR data within systematic uncertainties. 
For the case of a local source density "Ho = 10~ 5 Mpc -3 , 
the resultant ensemble fluctuations of the average mass 
composition on top of the ensemble mean were found to 
exist at a tolerable level for the analysis of present gen- 
eration CR composition data below 10 19 5 eV. 

However, the level of ensemble fluctuations for present 
generation spectral studies of the GZK suppression are 
potentially not ignorable if a smaller source density of 
Hq = 10~ 6 Mpc~ 3 is assumed, which still remains 
marginally consistent with angular correlation studies. 
In general, the amplitude of the ensemble fluctuations of 
the flux does not serve as a direct discriminator between 
CR models with different source compositions. This is 
apparent in the comparison between the iron and proton 
source models in Fig. [3] However, our study does pro- 
vide an analytic method to quantify these ensemble vari- 
ations enabling one to distinguish them from the GZK- 
suppression effect of the mean flux in a statistical way. 
Indeed, unless the actual source densities are much larger 
than those considered here, next generation experiments 
reaching accumulated exposures O(10 6 km 2 sryr) should 
be sufficiently sensitive to potentially discern these fluc- 
tuations. 

In addition to their spectral signatures, the hadronic 
interactions responsible for the GZK suppression are as- 
sociated with a high-energy flux of cosmogenic neutri- 
nos [55 63J and 7-rays [64-66 , following from photopion 
production of nucleons and subsequent decay. In addi- 
tion, photodisintcgration of high-energy nuclei can be fol- 
lowed by immediate photo-emission from excited daugh- 
ter nuclei |67j as well as neutrino-emission from /3-decay 
of the secondaries [HSJ HH] • 

While the presence of a suppression feature in the spec- 
trum is generally expected for all compositions, the flux 
of cosmogenic neutrinos and 7-rays subsequently pro- 
duced are both very sensitive to the CR source model. 
Indeed this difference permits information on these fluxes 
to be used as a probe of the composition or vice-versa. 
For instance, an upper limit on the proton fraction in 



UHE extra-galactic CRs can, in principle, be inferred 
from experimental bounds on both the diffuse flux of 
UHE neutrinos [35] and the diffuse flux of UHE pho- 
tons [23] . Furthermore, these two messengers starkly 
contrast in their subsequent propagation, with UHE neu- 
trinos freely propagating out to the Hubble-scale whereas 
UHE 7-rays being limited to tens of Mpc distance scales. 
This difference of scales highlights the fact that these two 
messengers can offer complementary information about 
the distant and local source distribution respectively. 

Moreover, the accompanying output into secondary 
electrons and positrons, in particular from Bethc-Hcitlcr 
pair production, feeds into electromagnetic cascades in 
the CMB and intergalactic magnetic fields. This leads 
to the accumulation of 7-rays at GeV-TeV energies. The 
observed diffuse 7-ray flux by Fermi-LAT |70j thus pro- 
vides a constraint on the total energy injected into such 
cascades over the Universe's entire history and can be 
translated into upper limits on the cosmic diffuse flux of 
photons and neutrinos [7iTf73] . 

In conclusion, future space-based observatories with 
colossal exposures, 0(1O 6 km 2 sryr), will provide the 
required large statistics at the high-energy end of the 
CR spectrum, allowing to separate ensemble fluctuation 
from the GZK suppression features on a statistical basis. 
In combination with information on the arrival-direction 
distribution of CRs and on the secondary fluxes of 7-rays 
and neutrinos these spectral features can provide a coher- 
ent picture for an indirect determination of the UHE CR 
nuclear composition [74] and will naturally complement 
the current direct measurements of EAS observables. 
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